#This code was used to conduct simulations of a stochastic game in the paper "Strategic Party Placement with a Dynamic Electorate."
#The game assumes that two parties place themselves on a single dimension and then compete in sequential elections with the same platform.
#This code allows three kinds of voter utilties: Quadratic (squared Euclidean distance), Absolute (Euclidean distance), and Directional.
#The actual simulations at the end of this file are all for Quadratic utilities.

#Code for running this document in R:
#source("/home/monogan/DISSERTATION/game/variance_SIMULATIONS_1000.R")
#source("F:/DISSERTATION/game/variance_SIMULATIONS_1000.R")
#source("Z:/variance_SIMULATIONS_1000.R")
#source("/netscr/monogan/variance_SIMULATIONS_1000.R")

#Set-up:
rm(list=ls())
options(width=1000)
library(xtable)

#Define Utility Functions:
Quadratic.A<-function(m.1,m.2,p,delta,theta.A,theta.D,shape){
    plogis(2*m.1*(theta.A-theta.D)+((theta.D)^2)-((theta.A)^2)+p)+delta*plogis(2*m.2*(theta.A-theta.D)+(theta.D^2)-(theta.A^2),s=shape)
    }
Quadratic.D<-function(m.1,m.2,p,delta,theta.A,theta.D,shape){
    (1-plogis(2*m.1*(theta.A-theta.D)+((theta.D)^2)-((theta.A)^2)+p))+delta*(1-plogis(2*m.2*(theta.A-theta.D)+(theta.D^2)-(theta.A^2),s=shape))
    }
    
Absolute.A<-function(m.1,m.2,v,delta,theta.A,theta.D,shape){
    plogis(-abs(m.1-theta.A)+abs(m.1-theta.D)+v)+delta*plogis(-abs(m.2-theta.A)+abs(m.2-theta.D),s=shape)
    }
Absolute.D<-function(m.1,m.2,v,delta,theta.A,theta.D,shape){
    (1-plogis(-abs(m.1-theta.A)+abs(m.1-theta.D)+v))+delta*(1-plogis(-abs(m.2-theta.A)+abs(m.2-theta.D),s=shape))
    }

Directional.A<-function(m.1,m.2,v,delta,theta.A,theta.D,shape){
    plogis(m.1*(theta.A-theta.D)+v)+delta*plogis(m.2*(theta.A-theta.D),s=shape)
    }
Directional.D<-function(m.1,m.2,v,delta,theta.A,theta.D,shape){
    (1-plogis(m.1*(theta.A-theta.D)+v))+delta*(1-plogis(m.2*(theta.A-theta.D),s=shape))
    }

# Define 'simulation' class
setClass('simulation', 
    representation(outcomeA='matrix', outcomeD='matrix', bestResponseA='numeric', bestResponseD='numeric', 
        equilibriumA='character', equilibriumD='character', parameters='character'))

# Simulation Function, which performs simulation and creates output
simulate<-function(type,v,delta,m.2,m.1=0.7,theta.A=seq(-1,1,.001),theta.D=seq(-1,1,.001),shape=1){
    #call the correct utility functions
    utilityA<-get(paste(type,"A",sep="."))
    utilityD<-get(paste(type,"D",sep="."))

    #create matrices and vectors
    outcomeA<-matrix(NA,2001,2001) #If NA isn't replaced, I need an error message that it's not numeric.
    outcomeD<-matrix(NA,2001,2001)
    bestResponseA<-rep(NA,2001) #Again, NA needs to be replaced.
    bestResponseD<-rep(NA,2001)
    equilibriumA<- 'NA' #NA may stick, but this is class 'character'
    equilibriumD<- 'NA'

    #matrix attributes
    rownames(outcomeA)<-seq(-1.0,1.0,0.001)
    colnames(outcomeA)<-seq(-1.0,1.0,0.001)
    rownames(outcomeD)<-seq(-1.0,1.0,0.001)
    colnames(outcomeD)<-seq(-1.0,1.0,0.001)
    names(bestResponseA)<-seq(-1.0,1.0,0.001)
    names(bestResponseD)<-seq(-1.0,1.0,0.001)

    #fill-in the utilities for all strategies for party A
    for (i in 1:2001){
            for (j in 1:2001){
                outcomeA[i,j]<-utilityA(m.1,m.2,v,delta,theta.A[i],theta.D[j],shape)
                }
        }

    #utilities for party D
    for (i in 1:2001){
            for (j in 1:2001){
                outcomeD[i,j]<-utilityD(m.1,m.2,v,delta,theta.A[i],theta.D[j],shape)
                }
        }

    #best responses for party A
    for (i in 1:2001){
        bestResponseA[i]<-which.max(outcomeA[,i])
        }

    #best responses for party D
    for (i in 1:2001){
        bestResponseD[i]<-which.max(outcomeD[i,])
        }

    #find the equilibria
    for (i in 1:2001){
        if (bestResponseD[bestResponseA[i]]==i){
            equilibriumA<-dimnames(outcomeA)[[1]][bestResponseA[i]]
            equilibriumD<-dimnames(outcomeD)[[2]][bestResponseD[bestResponseA[i]]]                  
                }
    }

    #save the input parameters in a slot that can be called
    parms<-c(type,v,delta,m.2,shape)
    
    #call the simulation class and fill the slots
    result<-new('simulation', outcomeA=outcomeA, outcomeD=outcomeD, bestResponseA=bestResponseA, bestResponseD=bestResponseD,
        equilibriumA=equilibriumA, equilibriumD=equilibriumD, parameters=parms)
    result
}

##################### RUN SIMULATIONS #########################################
#In these simulations, I have dropped movement of ideology and other voter calculuses from consideration.
#I simply manipulate the shape parameter

#Shape parameter is 1 by default. Variance=pi^2 / 3
treatment.1<-simulate("Quadratic",v=0.1,delta=0.0,m.2=-0.1)
treatment.2<-simulate("Quadratic",v=0.1,delta=0.1,m.2=-0.1)
treatment.3<-simulate("Quadratic",v=0.1,delta=0.5,m.2=-0.1)
treatment.4<-simulate("Quadratic",v=0.1,delta=0.9,m.2=-0.1)
treatment.5<-simulate("Quadratic",v=0.6,delta=0.0,m.2=-0.1)
treatment.6<-simulate("Quadratic",v=0.6,delta=0.1,m.2=-0.1)
treatment.7<-simulate("Quadratic",v=0.6,delta=0.5,m.2=-0.1)
treatment.8<-simulate("Quadratic",v=0.6,delta=0.9,m.2=-0.1)
treatment.9<-simulate("Quadratic",v=1.0,delta=0.0,m.2=-0.1)
treatment.10<-simulate("Quadratic",v=1.0,delta=0.1,m.2=-0.1)
treatment.11<-simulate("Quadratic",v=1.0,delta=0.5,m.2=-0.1)
treatment.12<-simulate("Quadratic",v=1.0,delta=0.9,m.2=-0.1)
treatment.13<-simulate("Quadratic",v=1.5,delta=0.0,m.2=-0.1)
treatment.14<-simulate("Quadratic",v=1.5,delta=0.1,m.2=-0.1)
treatment.15<-simulate("Quadratic",v=1.5,delta=0.5,m.2=-0.1)
treatment.16<-simulate("Quadratic",v=1.5,delta=0.9,m.2=-0.1)

#Shape parameter is 1/ sqrt(2). Variance=pi^2 / 6
treatment.17<-simulate("Quadratic",v=0.1,delta=0.0,m.2=-0.1,shape=1/sqrt(2))
treatment.18<-simulate("Quadratic",v=0.1,delta=0.1,m.2=-0.1,shape=1/sqrt(2))
treatment.19<-simulate("Quadratic",v=0.1,delta=0.5,m.2=-0.1,shape=1/sqrt(2))
treatment.20<-simulate("Quadratic",v=0.1,delta=0.9,m.2=-0.1,shape=1/sqrt(2))
treatment.21<-simulate("Quadratic",v=0.6,delta=0.0,m.2=-0.1,shape=1/sqrt(2))
treatment.22<-simulate("Quadratic",v=0.6,delta=0.1,m.2=-0.1,shape=1/sqrt(2))
treatment.23<-simulate("Quadratic",v=0.6,delta=0.5,m.2=-0.1,shape=1/sqrt(2))
treatment.24<-simulate("Quadratic",v=0.6,delta=0.9,m.2=-0.1,shape=1/sqrt(2))
treatment.25<-simulate("Quadratic",v=1.0,delta=0.0,m.2=-0.1,shape=1/sqrt(2))
treatment.26<-simulate("Quadratic",v=1.0,delta=0.1,m.2=-0.1,shape=1/sqrt(2))
treatment.27<-simulate("Quadratic",v=1.0,delta=0.5,m.2=-0.1,shape=1/sqrt(2))
treatment.28<-simulate("Quadratic",v=1.0,delta=0.9,m.2=-0.1,shape=1/sqrt(2))
treatment.29<-simulate("Quadratic",v=1.5,delta=0.0,m.2=-0.1,shape=1/sqrt(2))
treatment.30<-simulate("Quadratic",v=1.5,delta=0.1,m.2=-0.1,shape=1/sqrt(2))
treatment.31<-simulate("Quadratic",v=1.5,delta=0.5,m.2=-0.1,shape=1/sqrt(2))
treatment.32<-simulate("Quadratic",v=1.5,delta=0.9,m.2=-0.1,shape=1/sqrt(2))

#Shape parameter is 1/ sqrt(3). Variance=pi^2 / 9
treatment.33<-simulate("Quadratic",v=0.1,delta=0.0,m.2=-0.1,shape=1/sqrt(3))
treatment.34<-simulate("Quadratic",v=0.1,delta=0.1,m.2=-0.1,shape=1/sqrt(3))
treatment.35<-simulate("Quadratic",v=0.1,delta=0.5,m.2=-0.1,shape=1/sqrt(3))
treatment.36<-simulate("Quadratic",v=0.1,delta=0.9,m.2=-0.1,shape=1/sqrt(3))
treatment.37<-simulate("Quadratic",v=0.6,delta=0.0,m.2=-0.1,shape=1/sqrt(3))
treatment.38<-simulate("Quadratic",v=0.6,delta=0.1,m.2=-0.1,shape=1/sqrt(3))
treatment.39<-simulate("Quadratic",v=0.6,delta=0.5,m.2=-0.1,shape=1/sqrt(3))
treatment.40<-simulate("Quadratic",v=0.6,delta=0.9,m.2=-0.1,shape=1/sqrt(3))
treatment.41<-simulate("Quadratic",v=1.0,delta=0.0,m.2=-0.1,shape=1/sqrt(3))
treatment.42<-simulate("Quadratic",v=1.0,delta=0.1,m.2=-0.1,shape=1/sqrt(3))
treatment.43<-simulate("Quadratic",v=1.0,delta=0.5,m.2=-0.1,shape=1/sqrt(3))
treatment.44<-simulate("Quadratic",v=1.0,delta=0.9,m.2=-0.1,shape=1/sqrt(3))
treatment.45<-simulate("Quadratic",v=1.5,delta=0.0,m.2=-0.1,shape=1/sqrt(3))
treatment.46<-simulate("Quadratic",v=1.5,delta=0.1,m.2=-0.1,shape=1/sqrt(3))
treatment.47<-simulate("Quadratic",v=1.5,delta=0.5,m.2=-0.1,shape=1/sqrt(3))
treatment.48<-simulate("Quadratic",v=1.5,delta=0.9,m.2=-0.1,shape=1/sqrt(3))

##################### CREATE UTILITY TABLES #########################################
#Skipped for now

##################### PRINT BEST RESPONSES AND NASH EQUILIBRIA  ###########################
sink("/netscr/monogan/variance_equilibria_1000.txt")
#sink("Z:/variance_equilibria_1000.txt")
#sink("/home/monogan/DISSERTATION/game/variance_equilibria_1000.txt")
#sink("F:/DISSERTATION/game/variance_equilibria_1000.txt")

for(k in 1:48){

    treatment<-get(paste("treatment.",k,sep="")) #call the treatment

#    strategyA<-(treatment@bestResponseA-11)/10 #algebraically create strategies from indices
#    strategyD<-(treatment@bestResponseD-11)/10

    cat("TREATMENT ", k, ": \n", sep="")
    cat("utility=", treatment@parameters[1], ", v=", treatment@parameters[2], " , delta=", treatment@parameters[3], " , m.2=", treatment@parameters[4],
    " , shape=", treatment@parameters[5], "\n", sep="")
#    cat("Best Responses for Party A: \n")
#    print(strategyA)
#    cat("Best responses for Party D: \n")
#    print(strategyD)
    cat("Nash Equilibrium: ", treatment@equilibriumA, ",", treatment@equilibriumD, "\n \n")

}

sink()

